{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 求解一维线性波动方程"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. 定义单元算子和几何量"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "导入单元算子Operators和几何计算geomatry.jl，计算标准单元和全局单元上的节点分布、Vandermode矩阵和微分矩阵:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "using DifferentialEquations, Plots, BenchmarkTools, Sundials"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "include(\"geomatry.jl\")\n",
    "include(\"operators.jl\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "P = 8        # order of polynomials\n",
    "K = 10       # number of elements"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "9×9 Matrix{Float64}:\n",
       " 40.5      -1.84361    1.43139   -1.27425   …   1.43139   -1.84361    4.5\n",
       " -1.84361   6.79777   -0.586429   0.522048     -0.586429   0.755308  -1.84361\n",
       "  1.43139  -0.586429   4.09778   -0.405323      0.455309  -0.586429   1.43139\n",
       " -1.27425   0.522048  -0.405323   3.24742      -0.405323   0.522048  -1.27425\n",
       "  1.23047  -0.504111   0.391397  -0.348428      0.391397  -0.504111   1.23047\n",
       " -1.27425   0.522048  -0.405323   0.360825  …  -0.405323   0.522048  -1.27425\n",
       "  1.43139  -0.586429   0.455309  -0.405323      4.09778   -0.586429   1.43139\n",
       " -1.84361   0.755308  -0.586429   0.522048     -0.586429   6.79777   -1.84361\n",
       "  4.5      -1.84361    1.43139   -1.27425       1.43139   -1.84361   40.5"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "r = LegendreGaussLobatto(P)       # LGL solutions points at standard domain:[-1,1]\n",
    "h, x = Geomatry1D(0, 2π, P, K, r) # get element size h and global solution points x\n",
    "\n",
    "V = Vandermonde(P, r)      # Vandermonde matrix\n",
    "D = DiffMatrix(P, r, V)    # Differential operator\n",
    "Φ = V*V'                   # Lift operator"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. 初始化"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": "<?xml version=\"1.0\" encoding=\"utf-8\"?>\n<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"600\" height=\"400\" viewBox=\"0 0 2400 1600\">\n<defs>\n  <clipPath id=\"clip580\">\n    <rect x=\"0\" y=\"0\" width=\"2400\" height=\"1600\"/>\n  </clipPath>\n</defs>\n<path clip-path=\"url(#clip580)\" d=\"\nM0 1600 L2400 1600 L2400 0 L0 0  Z\n  \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n<defs>\n  <clipPath id=\"clip581\">\n    <rect x=\"480\" y=\"0\" width=\"1681\" height=\"1600\"/>\n  </clipPath>\n</defs>\n<path clip-path=\"url(#clip580)\" d=\"\nM193.936 1486.45 L2352.76 1486.45 L2352.76 47.2441 L193.936 47.2441  Z\n  \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n<defs>\n  <clipPath id=\"clip582\">\n    <rect x=\"193\" y=\"47\" width=\"2160\" height=\"1440\"/>\n  </clipPath>\n</defs>\n<polyline clip-path=\"url(#clip582)\" style=\"stroke:#000000; stroke-linecap:butt; stroke-linejoin:round; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n  255.035,1486.45 255.035,47.2441 \n  \"/>\n<polyline clip-path=\"url(#clip582)\" style=\"stroke:#000000; stroke-linecap:butt; stroke-linejoin:round; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n  579.173,1486.45 579.173,47.2441 \n  \"/>\n<polyline clip-path=\"url(#clip582)\" style=\"stroke:#000000; stroke-linecap:butt; stroke-linejoin:round; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n  903.312,1486.45 903.312,47.2441 \n  \"/>\n<polyline clip-path=\"url(#clip582)\" style=\"stroke:#000000; stroke-linecap:butt; stroke-linejoin:round; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n  1227.45,1486.45 1227.45,47.2441 \n  \"/>\n<polyline clip-path=\"url(#clip582)\" style=\"stroke:#000000; stroke-linecap:butt; stroke-linejoin:round; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n  1551.59,1486.45 1551.59,47.2441 \n  \"/>\n<polyline clip-path=\"url(#clip582)\" style=\"stroke:#000000; stroke-linecap:butt; stroke-linejoin:round; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n  1875.73,1486.45 1875.73,47.2441 \n  \"/>\n<polyline clip-path=\"url(#clip582)\" style=\"stroke:#000000; stroke-linecap:butt; stroke-linejoin:round; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n  2199.87,1486.45 2199.87,47.2441 \n  \"/>\n<polyline clip-path=\"url(#clip580)\" style=\"stroke:#000000; stroke-linecap:butt; stroke-linejoin:round; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n  193.936,1486.45 2352.76,1486.45 \n  \"/>\n<polyline clip-path=\"url(#clip580)\" style=\"stroke:#000000; stroke-linecap:butt; stroke-linejoin:round; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n  255.035,1486.45 255.035,1467.55 \n  \"/>\n<polyline clip-path=\"url(#clip580)\" style=\"stroke:#000000; stroke-linecap:butt; stroke-linejoin:round; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n  579.173,1486.45 579.173,1467.55 \n  \"/>\n<polyline clip-path=\"url(#clip580)\" style=\"stroke:#000000; stroke-linecap:butt; stroke-linejoin:round; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n  903.312,1486.45 903.312,1467.55 \n  \"/>\n<polyline clip-path=\"url(#clip580)\" style=\"stroke:#000000; stroke-linecap:butt; stroke-linejoin:round; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n  1227.45,1486.45 1227.45,1467.55 \n  \"/>\n<polyline clip-path=\"url(#clip580)\" style=\"stroke:#000000; stroke-linecap:butt; stroke-linejoin:round; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n  1551.59,1486.45 1551.59,1467.55 \n  \"/>\n<polyline clip-path=\"url(#clip580)\" style=\"stroke:#000000; stroke-linecap:butt; stroke-linejoin:round; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n  1875.73,1486.45 1875.73,1467.55 \n  \"/>\n<polyline clip-path=\"url(#clip580)\" style=\"stroke:#000000; stroke-linecap:butt; stroke-linejoin:round; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n  2199.87,1486.45 2199.87,1467.55 \n  \"/>\n<path clip-path=\"url(#clip580)\" d=\"M255.035 1515.64 Q251.424 1515.64 249.595 1519.2 Q247.789 1522.75 247.789 1529.87 Q247.789 1536.98 249.595 1540.55 Q251.424 1544.09 255.035 1544.09 Q258.669 1544.09 260.474 1540.55 Q262.303 1536.98 262.303 1529.87 Q262.303 1522.75 260.474 1519.2 Q258.669 1515.64 255.035 1515.64 M255.035 1511.93 Q260.845 1511.93 263.9 1516.54 Q266.979 1521.12 266.979 1529.87 Q266.979 1538.6 263.9 1543.21 Q260.845 1547.79 255.035 1547.79 Q249.225 1547.79 246.146 1543.21 Q243.09 1538.6 243.09 1529.87 Q243.09 1521.12 246.146 1516.54 Q249.225 1511.93 255.035 1511.93 Z\" fill=\"#000000\" fill-rule=\"evenodd\" fill-opacity=\"1\" /><path clip-path=\"url(#clip580)\" d=\"M569.555 1543.18 L577.194 1543.18 L577.194 1516.82 L568.884 1518.49 L568.884 1514.23 L577.148 1512.56 L581.824 1512.56 L581.824 1543.18 L589.463 1543.18 L589.463 1547.12 L569.555 1547.12 L569.555 1543.18 Z\" fill=\"#000000\" fill-rule=\"evenodd\" fill-opacity=\"1\" /><path clip-path=\"url(#clip580)\" d=\"M897.965 1543.18 L914.284 1543.18 L914.284 1547.12 L892.34 1547.12 L892.34 1543.18 Q895.002 1540.43 899.585 1535.8 Q904.191 1531.15 905.372 1529.81 Q907.617 1527.28 908.497 1525.55 Q909.4 1523.79 909.4 1522.1 Q909.4 1519.34 907.455 1517.61 Q905.534 1515.87 902.432 1515.87 Q900.233 1515.87 897.779 1516.63 Q895.349 1517.4 892.571 1518.95 L892.571 1514.23 Q895.395 1513.09 897.849 1512.51 Q900.303 1511.93 902.34 1511.93 Q907.71 1511.93 910.904 1514.62 Q914.099 1517.31 914.099 1521.8 Q914.099 1523.93 913.289 1525.85 Q912.502 1527.74 910.395 1530.34 Q909.816 1531.01 906.715 1534.23 Q903.613 1537.42 897.965 1543.18 Z\" fill=\"#000000\" fill-rule=\"evenodd\" fill-opacity=\"1\" /><path clip-path=\"url(#clip580)\" d=\"M1231.7 1528.49 Q1235.05 1529.2 1236.93 1531.47 Q1238.83 1533.74 1238.83 1537.07 Q1238.83 1542.19 1235.31 1544.99 Q1231.79 1547.79 1225.31 1547.79 Q1223.13 1547.79 1220.82 1547.35 Q1218.53 1546.93 1216.07 1546.08 L1216.07 1541.56 Q1218.02 1542.7 1220.33 1543.28 Q1222.65 1543.86 1225.17 1543.86 Q1229.57 1543.86 1231.86 1542.12 Q1234.17 1540.38 1234.17 1537.07 Q1234.17 1534.02 1232.02 1532.31 Q1229.89 1530.57 1226.07 1530.57 L1222.05 1530.57 L1222.05 1526.73 L1226.26 1526.73 Q1229.71 1526.73 1231.54 1525.36 Q1233.36 1523.97 1233.36 1521.38 Q1233.36 1518.72 1231.47 1517.31 Q1229.59 1515.87 1226.07 1515.87 Q1224.15 1515.87 1221.95 1516.29 Q1219.75 1516.7 1217.11 1517.58 L1217.11 1513.42 Q1219.78 1512.68 1222.09 1512.31 Q1224.43 1511.93 1226.49 1511.93 Q1231.81 1511.93 1234.92 1514.37 Q1238.02 1516.77 1238.02 1520.89 Q1238.02 1523.76 1236.37 1525.75 Q1234.73 1527.72 1231.7 1528.49 Z\" fill=\"#000000\" fill-rule=\"evenodd\" fill-opacity=\"1\" /><path clip-path=\"url(#clip580)\" d=\"M1554.6 1516.63 L1542.79 1535.08 L1554.6 1535.08 L1554.6 1516.63 M1553.37 1512.56 L1559.25 1512.56 L1559.25 1535.08 L1564.18 1535.08 L1564.18 1538.97 L1559.25 1538.97 L1559.25 1547.12 L1554.6 1547.12 L1554.6 1538.97 L1539 1538.97 L1539 1534.46 L1553.37 1512.56 Z\" fill=\"#000000\" fill-rule=\"evenodd\" fill-opacity=\"1\" /><path clip-path=\"url(#clip580)\" d=\"M1866.01 1512.56 L1884.36 1512.56 L1884.36 1516.5 L1870.29 1516.5 L1870.29 1524.97 Q1871.31 1524.62 1872.32 1524.46 Q1873.34 1524.27 1874.36 1524.27 Q1880.15 1524.27 1883.53 1527.44 Q1886.91 1530.62 1886.91 1536.03 Q1886.91 1541.61 1883.44 1544.71 Q1879.96 1547.79 1873.64 1547.79 Q1871.47 1547.79 1869.2 1547.42 Q1866.95 1547.05 1864.55 1546.31 L1864.55 1541.61 Q1866.63 1542.74 1868.85 1543.3 Q1871.07 1543.86 1873.55 1543.86 Q1877.56 1543.86 1879.89 1541.75 Q1882.23 1539.64 1882.23 1536.03 Q1882.23 1532.42 1879.89 1530.31 Q1877.56 1528.21 1873.55 1528.21 Q1871.68 1528.21 1869.8 1528.62 Q1867.95 1529.04 1866.01 1529.92 L1866.01 1512.56 Z\" fill=\"#000000\" fill-rule=\"evenodd\" fill-opacity=\"1\" /><path clip-path=\"url(#clip580)\" d=\"M2200.27 1527.98 Q2197.12 1527.98 2195.27 1530.13 Q2193.44 1532.28 2193.44 1536.03 Q2193.44 1539.76 2195.27 1541.93 Q2197.12 1544.09 2200.27 1544.09 Q2203.42 1544.09 2205.25 1541.93 Q2207.1 1539.76 2207.1 1536.03 Q2207.1 1532.28 2205.25 1530.13 Q2203.42 1527.98 2200.27 1527.98 M2209.55 1513.32 L2209.55 1517.58 Q2207.79 1516.75 2205.99 1516.31 Q2204.21 1515.87 2202.45 1515.87 Q2197.82 1515.87 2195.36 1519 Q2192.93 1522.12 2192.59 1528.44 Q2193.95 1526.43 2196.01 1525.36 Q2198.07 1524.27 2200.55 1524.27 Q2205.76 1524.27 2208.77 1527.44 Q2211.8 1530.59 2211.8 1536.03 Q2211.8 1541.36 2208.65 1544.57 Q2205.5 1547.79 2200.27 1547.79 Q2194.28 1547.79 2191.1 1543.21 Q2187.93 1538.6 2187.93 1529.87 Q2187.93 1521.68 2191.82 1516.82 Q2195.71 1511.93 2202.26 1511.93 Q2204.02 1511.93 2205.8 1512.28 Q2207.61 1512.63 2209.55 1513.32 Z\" fill=\"#000000\" fill-rule=\"evenodd\" fill-opacity=\"1\" /><polyline clip-path=\"url(#clip582)\" style=\"stroke:#000000; stroke-linecap:butt; stroke-linejoin:round; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n  193.936,1445.72 2352.76,1445.72 \n  \"/>\n<polyline clip-path=\"url(#clip582)\" style=\"stroke:#000000; stroke-linecap:butt; stroke-linejoin:round; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n  193.936,1106.28 2352.76,1106.28 \n  \"/>\n<polyline clip-path=\"url(#clip582)\" style=\"stroke:#000000; stroke-linecap:butt; stroke-linejoin:round; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n  193.936,766.846 2352.76,766.846 \n  \"/>\n<polyline clip-path=\"url(#clip582)\" style=\"stroke:#000000; stroke-linecap:butt; stroke-linejoin:round; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n  193.936,427.411 2352.76,427.411 \n  \"/>\n<polyline clip-path=\"url(#clip582)\" style=\"stroke:#000000; stroke-linecap:butt; stroke-linejoin:round; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n  193.936,87.9763 2352.76,87.9763 \n  \"/>\n<polyline clip-path=\"url(#clip580)\" style=\"stroke:#000000; stroke-linecap:butt; stroke-linejoin:round; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n  193.936,1486.45 193.936,47.2441 \n  \"/>\n<polyline clip-path=\"url(#clip580)\" style=\"stroke:#000000; stroke-linecap:butt; stroke-linejoin:round; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n  193.936,1445.72 212.834,1445.72 \n  \"/>\n<polyline clip-path=\"url(#clip580)\" style=\"stroke:#000000; stroke-linecap:butt; stroke-linejoin:round; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n  193.936,1106.28 212.834,1106.28 \n  \"/>\n<polyline clip-path=\"url(#clip580)\" style=\"stroke:#000000; stroke-linecap:butt; stroke-linejoin:round; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n  193.936,766.846 212.834,766.846 \n  \"/>\n<polyline clip-path=\"url(#clip580)\" style=\"stroke:#000000; stroke-linecap:butt; stroke-linejoin:round; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n  193.936,427.411 212.834,427.411 \n  \"/>\n<polyline clip-path=\"url(#clip580)\" style=\"stroke:#000000; stroke-linecap:butt; stroke-linejoin:round; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n  193.936,87.9763 212.834,87.9763 \n  \"/>\n<path clip-path=\"url(#clip580)\" d=\"M50.9921 1446.17 L80.6679 1446.17 L80.6679 1450.1 L50.9921 1450.1 L50.9921 1446.17 Z\" fill=\"#000000\" fill-rule=\"evenodd\" fill-opacity=\"1\" /><path clip-path=\"url(#clip580)\" d=\"M91.5706 1459.06 L99.2095 1459.06 L99.2095 1432.69 L90.8993 1434.36 L90.8993 1430.1 L99.1632 1428.44 L103.839 1428.44 L103.839 1459.06 L111.478 1459.06 L111.478 1463 L91.5706 1463 L91.5706 1459.06 Z\" fill=\"#000000\" fill-rule=\"evenodd\" fill-opacity=\"1\" /><path clip-path=\"url(#clip580)\" d=\"M120.922 1457.12 L125.807 1457.12 L125.807 1463 L120.922 1463 L120.922 1457.12 Z\" fill=\"#000000\" fill-rule=\"evenodd\" fill-opacity=\"1\" /><path clip-path=\"url(#clip580)\" d=\"M145.992 1431.51 Q142.381 1431.51 140.552 1435.08 Q138.746 1438.62 138.746 1445.75 Q138.746 1452.86 140.552 1456.42 Q142.381 1459.96 145.992 1459.96 Q149.626 1459.96 151.431 1456.42 Q153.26 1452.86 153.26 1445.75 Q153.26 1438.62 151.431 1435.08 Q149.626 1431.51 145.992 1431.51 M145.992 1427.81 Q151.802 1427.81 154.857 1432.42 Q157.936 1437 157.936 1445.75 Q157.936 1454.48 154.857 1459.08 Q151.802 1463.67 145.992 1463.67 Q140.181 1463.67 137.103 1459.08 Q134.047 1454.48 134.047 1445.75 Q134.047 1437 137.103 1432.42 Q140.181 1427.81 145.992 1427.81 Z\" fill=\"#000000\" fill-rule=\"evenodd\" fill-opacity=\"1\" /><path clip-path=\"url(#clip580)\" d=\"M51.9875 1106.73 L81.6633 1106.73 L81.6633 1110.67 L51.9875 1110.67 L51.9875 1106.73 Z\" fill=\"#000000\" fill-rule=\"evenodd\" fill-opacity=\"1\" /><path clip-path=\"url(#clip580)\" d=\"M101.756 1092.08 Q98.1447 1092.08 96.316 1095.64 Q94.5104 1099.19 94.5104 1106.32 Q94.5104 1113.42 96.316 1116.99 Q98.1447 1120.53 101.756 1120.53 Q105.39 1120.53 107.196 1116.99 Q109.024 1113.42 109.024 1106.32 Q109.024 1099.19 107.196 1095.64 Q105.39 1092.08 101.756 1092.08 M101.756 1088.38 Q107.566 1088.38 110.621 1092.98 Q113.7 1097.57 113.7 1106.32 Q113.7 1115.04 110.621 1119.65 Q107.566 1124.23 101.756 1124.23 Q95.9456 1124.23 92.8669 1119.65 Q89.8114 1115.04 89.8114 1106.32 Q89.8114 1097.57 92.8669 1092.98 Q95.9456 1088.38 101.756 1088.38 Z\" fill=\"#000000\" fill-rule=\"evenodd\" fill-opacity=\"1\" /><path clip-path=\"url(#clip580)\" d=\"M121.918 1117.68 L126.802 1117.68 L126.802 1123.56 L121.918 1123.56 L121.918 1117.68 Z\" fill=\"#000000\" fill-rule=\"evenodd\" fill-opacity=\"1\" /><path clip-path=\"url(#clip580)\" d=\"M137.033 1089 L155.39 1089 L155.39 1092.94 L141.316 1092.94 L141.316 1101.41 Q142.334 1101.06 143.353 1100.9 Q144.371 1100.71 145.39 1100.71 Q151.177 1100.71 154.556 1103.89 Q157.936 1107.06 157.936 1112.47 Q157.936 1118.05 154.464 1121.15 Q150.992 1124.23 144.672 1124.23 Q142.496 1124.23 140.228 1123.86 Q137.982 1123.49 135.575 1122.75 L135.575 1118.05 Q137.658 1119.19 139.881 1119.74 Q142.103 1120.3 144.58 1120.3 Q148.584 1120.3 150.922 1118.19 Q153.26 1116.08 153.26 1112.47 Q153.26 1108.86 150.922 1106.76 Q148.584 1104.65 144.58 1104.65 Q142.705 1104.65 140.83 1105.07 Q138.978 1105.48 137.033 1106.36 L137.033 1089 Z\" fill=\"#000000\" fill-rule=\"evenodd\" fill-opacity=\"1\" /><path clip-path=\"url(#clip580)\" d=\"M100.76 752.645 Q97.1493 752.645 95.3206 756.209 Q93.515 759.751 93.515 766.881 Q93.515 773.987 95.3206 777.552 Q97.1493 781.094 100.76 781.094 Q104.395 781.094 106.2 777.552 Q108.029 773.987 108.029 766.881 Q108.029 759.751 106.2 756.209 Q104.395 752.645 100.76 752.645 M100.76 748.941 Q106.571 748.941 109.626 753.547 Q112.705 758.131 112.705 766.881 Q112.705 775.608 109.626 780.214 Q106.571 784.797 100.76 784.797 Q94.9502 784.797 91.8715 780.214 Q88.816 775.608 88.816 766.881 Q88.816 758.131 91.8715 753.547 Q94.9502 748.941 100.76 748.941 Z\" fill=\"#000000\" fill-rule=\"evenodd\" fill-opacity=\"1\" /><path clip-path=\"url(#clip580)\" d=\"M120.922 778.246 L125.807 778.246 L125.807 784.126 L120.922 784.126 L120.922 778.246 Z\" fill=\"#000000\" fill-rule=\"evenodd\" fill-opacity=\"1\" /><path clip-path=\"url(#clip580)\" d=\"M145.992 752.645 Q142.381 752.645 140.552 756.209 Q138.746 759.751 138.746 766.881 Q138.746 773.987 140.552 777.552 Q142.381 781.094 145.992 781.094 Q149.626 781.094 151.431 777.552 Q153.26 773.987 153.26 766.881 Q153.26 759.751 151.431 756.209 Q149.626 752.645 145.992 752.645 M145.992 748.941 Q151.802 748.941 154.857 753.547 Q157.936 758.131 157.936 766.881 Q157.936 775.608 154.857 780.214 Q151.802 784.797 145.992 784.797 Q140.181 784.797 137.103 780.214 Q134.047 775.608 134.047 766.881 Q134.047 758.131 137.103 753.547 Q140.181 748.941 145.992 748.941 Z\" fill=\"#000000\" fill-rule=\"evenodd\" fill-opacity=\"1\" /><path clip-path=\"url(#clip580)\" d=\"M101.756 413.21 Q98.1447 413.21 96.316 416.775 Q94.5104 420.316 94.5104 427.446 Q94.5104 434.552 96.316 438.117 Q98.1447 441.659 101.756 441.659 Q105.39 441.659 107.196 438.117 Q109.024 434.552 109.024 427.446 Q109.024 420.316 107.196 416.775 Q105.39 413.21 101.756 413.21 M101.756 409.506 Q107.566 409.506 110.621 414.113 Q113.7 418.696 113.7 427.446 Q113.7 436.173 110.621 440.779 Q107.566 445.362 101.756 445.362 Q95.9456 445.362 92.8669 440.779 Q89.8114 436.173 89.8114 427.446 Q89.8114 418.696 92.8669 414.113 Q95.9456 409.506 101.756 409.506 Z\" fill=\"#000000\" fill-rule=\"evenodd\" fill-opacity=\"1\" /><path clip-path=\"url(#clip580)\" d=\"M121.918 438.812 L126.802 438.812 L126.802 444.691 L121.918 444.691 L121.918 438.812 Z\" fill=\"#000000\" fill-rule=\"evenodd\" fill-opacity=\"1\" /><path clip-path=\"url(#clip580)\" d=\"M137.033 410.131 L155.39 410.131 L155.39 414.066 L141.316 414.066 L141.316 422.538 Q142.334 422.191 143.353 422.029 Q144.371 421.844 145.39 421.844 Q151.177 421.844 154.556 425.015 Q157.936 428.187 157.936 433.603 Q157.936 439.182 154.464 442.284 Q150.992 445.362 144.672 445.362 Q142.496 445.362 140.228 444.992 Q137.982 444.622 135.575 443.881 L135.575 439.182 Q137.658 440.316 139.881 440.872 Q142.103 441.427 144.58 441.427 Q148.584 441.427 150.922 439.321 Q153.26 437.214 153.26 433.603 Q153.26 429.992 150.922 427.886 Q148.584 425.779 144.58 425.779 Q142.705 425.779 140.83 426.196 Q138.978 426.613 137.033 427.492 L137.033 410.131 Z\" fill=\"#000000\" fill-rule=\"evenodd\" fill-opacity=\"1\" /><path clip-path=\"url(#clip580)\" d=\"M91.5706 101.321 L99.2095 101.321 L99.2095 74.9555 L90.8993 76.6222 L90.8993 72.3629 L99.1632 70.6963 L103.839 70.6963 L103.839 101.321 L111.478 101.321 L111.478 105.256 L91.5706 105.256 L91.5706 101.321 Z\" fill=\"#000000\" fill-rule=\"evenodd\" fill-opacity=\"1\" /><path clip-path=\"url(#clip580)\" d=\"M120.922 99.3767 L125.807 99.3767 L125.807 105.256 L120.922 105.256 L120.922 99.3767 Z\" fill=\"#000000\" fill-rule=\"evenodd\" fill-opacity=\"1\" /><path clip-path=\"url(#clip580)\" d=\"M145.992 73.775 Q142.381 73.775 140.552 77.3398 Q138.746 80.8814 138.746 88.011 Q138.746 95.1174 140.552 98.6822 Q142.381 102.224 145.992 102.224 Q149.626 102.224 151.431 98.6822 Q153.26 95.1174 153.26 88.011 Q153.26 80.8814 151.431 77.3398 Q149.626 73.775 145.992 73.775 M145.992 70.0713 Q151.802 70.0713 154.857 74.6777 Q157.936 79.261 157.936 88.011 Q157.936 96.7378 154.857 101.344 Q151.802 105.928 145.992 105.928 Q140.181 105.928 137.103 101.344 Q134.047 96.7378 134.047 88.011 Q134.047 79.261 137.103 74.6777 Q140.181 70.0713 145.992 70.0713 Z\" fill=\"#000000\" fill-rule=\"evenodd\" fill-opacity=\"1\" /><polyline clip-path=\"url(#clip582)\" style=\"stroke:#009af9; stroke-linecap:butt; stroke-linejoin:round; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n  255.035,766.846 265.242,745.471 287.907,698.116 319.889,631.92 356.866,557.064 393.842,484.934 425.824,425.47 448.489,385.307 458.697,367.816 \n  468.905,350.721 491.569,314.263 523.551,266.62 560.528,217.629 597.505,175.777 629.487,145.759 652.151,128.128 662.359,121.203 672.567,114.917 \n  695.232,103.281 727.214,92.3887 764.19,87.9763 801.167,92.3887 833.149,103.281 855.814,114.917 866.021,121.203 876.229,128.128 898.894,145.759 \n  930.876,175.777 967.853,217.629 1004.83,266.62 1036.81,314.263 1059.48,350.721 1069.68,367.816 1079.89,385.307 1102.56,425.47 1134.54,484.934 \n  1171.51,557.064 1208.49,631.92 1240.47,698.116 1263.14,745.471 1273.35,766.846 1283.55,788.221 1306.22,835.576 1338.2,901.771 1375.18,976.628 \n  1412.15,1048.76 1444.14,1108.22 1466.8,1148.38 1477.01,1165.88 1487.22,1182.97 1509.88,1219.43 1541.86,1267.07 1578.84,1316.06 1615.82,1357.91 \n  1647.8,1387.93 1670.46,1405.56 1680.67,1412.49 1690.88,1418.77 1713.54,1430.41 1745.52,1441.3 1782.5,1445.72 1819.48,1441.3 1851.46,1430.41 \n  1874.12,1418.77 1884.33,1412.49 1894.54,1405.56 1917.21,1387.93 1949.19,1357.91 1986.16,1316.06 2023.14,1267.07 2055.12,1219.43 2077.79,1182.97 \n  2087.99,1165.88 2098.2,1148.38 2120.87,1108.22 2152.85,1048.76 2189.83,976.628 2226.8,901.771 2258.78,835.576 2281.45,788.221 2291.66,766.846 \n  \n  \"/>\n<path clip-path=\"url(#clip580)\" d=\"\nM1914.61 198.898 L2280.8 198.898 L2280.8 95.2176 L1914.61 95.2176  Z\n  \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n<polyline clip-path=\"url(#clip580)\" style=\"stroke:#000000; stroke-linecap:butt; stroke-linejoin:round; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n  1914.61,198.898 2280.8,198.898 2280.8,95.2176 1914.61,95.2176 1914.61,198.898 \n  \"/>\n<polyline clip-path=\"url(#clip580)\" style=\"stroke:#009af9; stroke-linecap:butt; stroke-linejoin:round; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n  1938.59,147.058 2082.52,147.058 \n  \"/>\n<path clip-path=\"url(#clip580)\" d=\"M2124.93 139.176 L2124.93 143.203 Q2123.12 142.277 2121.18 141.815 Q2119.23 141.352 2117.15 141.352 Q2113.98 141.352 2112.38 142.324 Q2110.81 143.296 2110.81 145.24 Q2110.81 146.722 2111.94 147.578 Q2113.08 148.412 2116.5 149.176 L2117.96 149.5 Q2122.5 150.472 2124.4 152.254 Q2126.32 154.014 2126.32 157.185 Q2126.32 160.796 2123.45 162.902 Q2120.6 165.009 2115.6 165.009 Q2113.52 165.009 2111.25 164.592 Q2109 164.199 2106.5 163.388 L2106.5 158.99 Q2108.86 160.217 2111.16 160.842 Q2113.45 161.444 2115.69 161.444 Q2118.7 161.444 2120.32 160.426 Q2121.94 159.384 2121.94 157.509 Q2121.94 155.773 2120.76 154.847 Q2119.6 153.921 2115.65 153.064 L2114.16 152.717 Q2110.21 151.884 2108.45 150.171 Q2106.69 148.435 2106.69 145.426 Q2106.69 141.768 2109.28 139.778 Q2111.87 137.787 2116.64 137.787 Q2119 137.787 2121.09 138.134 Q2123.17 138.481 2124.93 139.176 Z\" fill=\"#000000\" fill-rule=\"evenodd\" fill-opacity=\"1\" /><path clip-path=\"url(#clip580)\" d=\"M2133.1 138.412 L2137.36 138.412 L2137.36 164.338 L2133.1 164.338 L2133.1 138.412 M2133.1 128.319 L2137.36 128.319 L2137.36 133.713 L2133.1 133.713 L2133.1 128.319 Z\" fill=\"#000000\" fill-rule=\"evenodd\" fill-opacity=\"1\" /><path clip-path=\"url(#clip580)\" d=\"M2167.82 148.689 L2167.82 164.338 L2163.56 164.338 L2163.56 148.828 Q2163.56 145.148 2162.13 143.319 Q2160.69 141.49 2157.82 141.49 Q2154.37 141.49 2152.38 143.69 Q2150.39 145.889 2150.39 149.685 L2150.39 164.338 L2146.11 164.338 L2146.11 138.412 L2150.39 138.412 L2150.39 142.44 Q2151.92 140.102 2153.98 138.944 Q2156.06 137.787 2158.77 137.787 Q2163.24 137.787 2165.53 140.565 Q2167.82 143.319 2167.82 148.689 Z\" fill=\"#000000\" fill-rule=\"evenodd\" fill-opacity=\"1\" /><path clip-path=\"url(#clip580)\" d=\"M2186.55 128.366 Q2183.45 133.69 2181.94 138.898 Q2180.44 144.106 2180.44 149.453 Q2180.44 154.801 2181.94 160.055 Q2183.47 165.287 2186.55 170.588 L2182.84 170.588 Q2179.37 165.148 2177.64 159.893 Q2175.92 154.639 2175.92 149.453 Q2175.92 144.291 2177.64 139.06 Q2179.35 133.828 2182.84 128.366 L2186.55 128.366 Z\" fill=\"#000000\" fill-rule=\"evenodd\" fill-opacity=\"1\" /><path clip-path=\"url(#clip580)\" d=\"M2216.36 138.412 L2206.99 151.027 L2216.85 164.338 L2211.83 164.338 L2204.28 154.152 L2196.73 164.338 L2191.71 164.338 L2201.78 150.773 L2192.57 138.412 L2197.59 138.412 L2204.47 147.648 L2211.34 138.412 L2216.36 138.412 Z\" fill=\"#000000\" fill-rule=\"evenodd\" fill-opacity=\"1\" /><path clip-path=\"url(#clip580)\" d=\"M2222.2 128.366 L2225.9 128.366 Q2229.37 133.828 2231.09 139.06 Q2232.82 144.291 2232.82 149.453 Q2232.82 154.639 2231.09 159.893 Q2229.37 165.148 2225.9 170.588 L2222.2 170.588 Q2225.28 165.287 2226.78 160.055 Q2228.31 154.801 2228.31 149.453 Q2228.31 144.106 2226.78 138.898 Q2225.28 133.69 2222.2 128.366 Z\" fill=\"#000000\" fill-rule=\"evenodd\" fill-opacity=\"1\" /></svg>\n"
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@inbounds u0 = sin.(x) \n",
    "plot(x[:], u0[:],label=\"sin(x)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. 右端项函数"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.1 基础版"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "rhs (generic function with 1 method)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function rhs(u,p,t)\n",
    "    D, Φ, h = p\n",
    "    Np, K = size(u)\n",
    "    α = 0\n",
    "    du = zeros(Np, K)\n",
    "    du[1,1] = (α/2.0 - 1)*(u[1,1]-u[end,end])\n",
    "    du[1,2:K] = (α/2.0 - 1)*(u[1,2:K]-u[end,1:K-1])\n",
    "    du[end,end] = α/2.0*(u[end,end]-u[1,1])\n",
    "    du[end,1:K-1]=α/2.0*(u[end,1:K-1]-u[1,2:K])\n",
    "    \n",
    "    return -2.0/h*(D*u .- Φ*du) \n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- 显式时间格式"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\u001b[36mODEProblem\u001b[0m with uType \u001b[36mMatrix{Float64}\u001b[0m and tType \u001b[36mFloat64\u001b[0m. In-place: \u001b[36mfalse\u001b[0m\n",
       "timespan: (0.0, 1.0)\n",
       "u0: 9×10 Matrix{Float64}:\n",
       " 0.0        0.587785  0.951057  …  -0.951057  -0.951057  -0.587785\n",
       " 0.0314867  0.612967  0.960315     -0.960315  -0.940855  -0.56202\n",
       " 0.101241   0.666671  0.977455     -0.977455  -0.914885  -0.502859\n",
       " 0.19875    0.736851  0.9935       -0.9935    -0.870666  -0.415267\n",
       " 0.309017   0.809017  1.0          -1.0       -0.809017  -0.309017\n",
       " 0.415267   0.870666  0.9935    …  -0.9935    -0.736851  -0.19875\n",
       " 0.502859   0.914885  0.977455     -0.977455  -0.666671  -0.101241\n",
       " 0.56202    0.940855  0.960315     -0.960315  -0.612967  -0.0314867\n",
       " 0.587785   0.951057  0.951057     -0.951057  -0.587785  -2.44929e-16"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "prob = ODEProblem(rhs, u0, (0.0,1.0),(D, Φ, h))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  1.967589 seconds (5.63 M allocations: 326.794 MiB, 2.69% gc time, 99.97% compilation time)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "retcode: Success\n",
       "Interpolation: 1st order linear\n",
       "t: 2-element Vector{Float64}:\n",
       " 0.0\n",
       " 1.0\n",
       "u: 2-element Vector{Matrix{Float64}}:\n",
       " [0.0 0.5877852522924731 … -0.9510565162951536 -0.5877852522924734; 0.031486749455803056 0.6129671257739798 … -0.9408550126738896 -0.5620204949592386; … ; 0.5620204949592384 0.9408550126738895 … -0.6129671257739798 -0.03148674945580308; 0.5877852522924731 0.9510565162951535 … -0.5877852522924734 -2.4492935982947064e-16]\n",
       " [-0.8414018239068098 -0.36325380133230345 … -0.7738178938452396 -0.9984143611514931; -0.8240522731953155 -0.33365452928529354 … -0.793455206538053 -0.9996505294823455; … ; -0.3923424140010108 0.22325101024290656 … -0.996044018680173 -0.8580630168071015; -0.3631466346382578 0.2537920238242879 … -0.998310871780859 -0.8415061948317838]"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@time u1 = solve(prob,Tsit5(),save_everystep=false)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- BDF隐式时间格式"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  1.419134 seconds (4.55 M allocations: 286.249 MiB, 3.49% gc time)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "retcode: Success\n",
       "Interpolation: 1st order linear\n",
       "t: 2-element Vector{Float64}:\n",
       " 0.0\n",
       " 1.0\n",
       "u: 2-element Vector{Matrix{Float64}}:\n",
       " [0.0 0.5877852522924731 … -0.9510565162951536 -0.5877852522924734; 0.031486749455803056 0.6129671257739798 … -0.9408550126738896 -0.5620204949592386; … ; 0.5620204949592384 0.9408550126738895 … -0.6129671257739798 -0.03148674945580308; 0.5877852522924731 0.9510565162951535 … -0.5877852522924734 -2.4492935982947064e-16]\n",
       " [-0.8414191423804401 -0.3631065856821248 … -0.7739486342866062 -0.9983545188792516; -0.823992142791415 -0.3335934892795848 … -0.7934784812904109 -0.9996510670860366; … ; -0.3922600327893867 0.22332530739193066 … -0.9960361868050959 -0.8580143066961294; -0.3631052455724007 0.25390042410751157 … -0.9983422678399212 -0.8414241644701116]"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@time solve(prob,CVODE_BDF(linear_solver=:GMRES),save_everystep=false)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- TRBDF2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  6.157647 seconds (23.02 M allocations: 1.455 GiB, 4.93% gc time, 99.50% compilation time)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "retcode: Success\n",
       "Interpolation: 1st order linear\n",
       "t: 1-element Vector{Float64}:\n",
       " 1.0\n",
       "u: 1-element Vector{Matrix{Float64}}:\n",
       " [-0.840206829147325 -0.3610671326394734 … -0.7752643143090131 -0.9984160745006739; -0.8227193289432748 -0.3315273730099106 … -0.7947652797007239 -0.999660464421813; … ; -0.39024883599301596 0.22542524810245224 … -0.996181594058928 -0.8568611288092353; -0.36106713264120766 0.2559879363147554 … -0.9984160745011579 -0.8402068291486959]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@time sol = solve(prob,TRBDF2(autodiff=false),save_everystep=false,save_start=false)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Rosenbrock23"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  2.688094 seconds (8.07 M allocations: 569.994 MiB, 2.90% gc time, 99.28% compilation time)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "retcode: Success\n",
       "Interpolation: 1st order linear\n",
       "t: 1-element Vector{Float64}:\n",
       " 1.0\n",
       "u: 1-element Vector{Matrix{Float64}}:\n",
       " [-0.8411880433553592 -0.36270119332706596 … -0.7742079851434087 -0.9983696512225814; -0.8237450669454046 -0.3331788780363053 … -0.7937527672206862 -0.9996686379158554; … ; -0.3918638323570215 0.2237478461817207 … -0.9960806200032617 -0.8577968452699789; -0.36270119322441696 0.25432518450636127 … -0.9983696512048045 -0.8411880432675355]"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@time sol = solve(prob,Rosenbrock23(autodiff=false),save_everystep=false,save_start=false)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.2 Update1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "rhs_update1 (generic function with 1 method)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function rhs_update1(u,p,t)\n",
    "   D, Φ, h, df, Np, K, α = p\n",
    "    df[1,1] = (α/2.0 - 1)*(u[1,1]-u[end,end])\n",
    "    df[1,2:K] = (α/2.0 - 1)*(u[1,2:K]-u[end,1:K-1])\n",
    "    df[end,end] = α/2.0*(u[end,end]-u[1,1])\n",
    "    df[end,1:K-1]=α/2.0*(u[end,1:K-1]-u[1,2:K])\n",
    "    return -2.0/h*(D*u .- Φ*df)  \n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 显式时间格式"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  1.028176 seconds (2.34 M allocations: 137.787 MiB, 2.61% gc time, 99.95% compilation time)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "retcode: Success\n",
       "Interpolation: 1st order linear\n",
       "t: 2-element Vector{Float64}:\n",
       " 0.0\n",
       " 1.0\n",
       "u: 2-element Vector{Matrix{Float64}}:\n",
       " [0.0 0.5877852522924731 … -0.9510565162951536 -0.5877852522924734; 0.031486749455803056 0.6129671257739798 … -0.9408550126738896 -0.5620204949592386; … ; 0.5620204949592384 0.9408550126738895 … -0.6129671257739798 -0.03148674945580308; 0.5877852522924731 0.9510565162951535 … -0.5877852522924734 -2.4492935982947064e-16]\n",
       " [-0.8414018239068098 -0.36325380133230345 … -0.7738178938452396 -0.9984143611514931; -0.8240522731953155 -0.33365452928529354 … -0.793455206538053 -0.9996505294823455; … ; -0.3923424140010108 0.22325101024290656 … -0.996044018680173 -0.8580630168071015; -0.3631466346382578 0.2537920238242879 … -0.998310871780859 -0.8415061948317838]"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Np = P+1\n",
    "df = zeros(Np, K)\n",
    "p = (D, Φ, h, df, Np, K, 0)\n",
    "prob = ODEProblem(rhs_update1, u0, (0.0,1.0), p)\n",
    "@time solve(prob,Tsit5(),save_everystep=false)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 隐式时间格式"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  0.524892 seconds (1.57 M allocations: 100.434 MiB, 4.51% gc time, 99.82% compilation time)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "retcode: Success\n",
       "Interpolation: 1st order linear\n",
       "t: 2-element Vector{Float64}:\n",
       " 0.0\n",
       " 1.0\n",
       "u: 2-element Vector{Matrix{Float64}}:\n",
       " [0.0 0.5877852522924731 … -0.9510565162951536 -0.5877852522924734; 0.031486749455803056 0.6129671257739798 … -0.9408550126738896 -0.5620204949592386; … ; 0.5620204949592384 0.9408550126738895 … -0.6129671257739798 -0.03148674945580308; 0.5877852522924731 0.9510565162951535 … -0.5877852522924734 -2.4492935982947064e-16]\n",
       " [-0.8414191423804401 -0.3631065856821248 … -0.7739486342866062 -0.9983545188792516; -0.823992142791415 -0.3335934892795848 … -0.7934784812904109 -0.9996510670860366; … ; -0.3922600327893867 0.22332530739193066 … -0.9960361868050959 -0.8580143066961294; -0.3631052455724007 0.25390042410751157 … -0.9983422678399212 -0.8414241644701116]"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Np = P+1\n",
    "df = zeros(Np, K)\n",
    "p = (D, Φ, h, df, Np, K, 0)\n",
    "prob = ODEProblem(rhs_update1, u0, (0.0,1.0), p)\n",
    "@time solve(prob,CVODE_BDF(linear_solver=:GMRES),save_everystep=false)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "BenchmarkTools.Trial: 10000 samples with 1 evaluation.\n",
       " Range \u001b[90m(\u001b[39m\u001b[36m\u001b[1mmin\u001b[22m\u001b[39m … \u001b[35mmax\u001b[39m\u001b[90m):  \u001b[39m\u001b[36m\u001b[1m219.561 μs\u001b[22m\u001b[39m … \u001b[35m 13.006 ms\u001b[39m  \u001b[90m┊\u001b[39m GC \u001b[90m(\u001b[39mmin … max\u001b[90m): \u001b[39m 0.00% … 80.51%\n",
       " Time  \u001b[90m(\u001b[39m\u001b[34m\u001b[1mmedian\u001b[22m\u001b[39m\u001b[90m):     \u001b[39m\u001b[34m\u001b[1m239.583 μs               \u001b[22m\u001b[39m\u001b[90m┊\u001b[39m GC \u001b[90m(\u001b[39mmedian\u001b[90m):    \u001b[39m 0.00%\n",
       " Time  \u001b[90m(\u001b[39m\u001b[32m\u001b[1mmean\u001b[22m\u001b[39m ± \u001b[32mσ\u001b[39m\u001b[90m):   \u001b[39m\u001b[32m\u001b[1m288.508 μs\u001b[22m\u001b[39m ± \u001b[32m769.181 μs\u001b[39m  \u001b[90m┊\u001b[39m GC \u001b[90m(\u001b[39mmean ± σ\u001b[90m):  \u001b[39m13.66% ±  5.01%\n",
       "\n",
       "  \u001b[39m▂\u001b[39m▆\u001b[39m█\u001b[39m▇\u001b[39m▆\u001b[39m▄\u001b[39m▅\u001b[34m▇\u001b[39m\u001b[39m█\u001b[39m▇\u001b[39m▆\u001b[39m▃\u001b[39m▁\u001b[39m▁\u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[32m \u001b[39m\u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m▂\n",
       "  \u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[34m█\u001b[39m\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m▇\u001b[39m█\u001b[39m█\u001b[39m▇\u001b[39m█\u001b[39m▇\u001b[32m▇\u001b[39m\u001b[39m▇\u001b[39m▇\u001b[39m▆\u001b[39m▆\u001b[39m▆\u001b[39m▆\u001b[39m▅\u001b[39m▆\u001b[39m▆\u001b[39m▅\u001b[39m▆\u001b[39m▅\u001b[39m▅\u001b[39m▅\u001b[39m▅\u001b[39m▅\u001b[39m▄\u001b[39m▅\u001b[39m▄\u001b[39m▄\u001b[39m▄\u001b[39m▄\u001b[39m▄\u001b[39m▄\u001b[39m▄\u001b[39m▂\u001b[39m▅\u001b[39m▃\u001b[39m▃\u001b[39m▅\u001b[39m▃\u001b[39m▂\u001b[39m▃\u001b[39m▅\u001b[39m▅\u001b[39m \u001b[39m█\n",
       "  220 μs\u001b[90m        \u001b[39m\u001b[90mHistogram: \u001b[39m\u001b[90m\u001b[1mlog(\u001b[22m\u001b[39m\u001b[90mfrequency\u001b[39m\u001b[90m\u001b[1m)\u001b[22m\u001b[39m\u001b[90m by time\u001b[39m        386 μs \u001b[0m\u001b[1m<\u001b[22m\n",
       "\n",
       " Memory estimate\u001b[90m: \u001b[39m\u001b[33m392.19 KiB\u001b[39m, allocs estimate\u001b[90m: \u001b[39m\u001b[33m1611\u001b[39m."
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@benchmark solve(prob,CVODE_BDF(linear_solver=:GMRES),save_everystep=false)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Update3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "rhs_update3 (generic function with 1 method)"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function rhs_update3(u,p,t)\n",
    "    D, Φ, h, df, Np, dF_d, dF_c, K, α = p\n",
    "    df[1,1] = (α/2.0 - 1)*(u[1,1]-u[end,end])\n",
    "    df[1,2:K] = (α/2.0 - 1)*(u[1,2:K]-u[end,1:K-1])\n",
    "    df[end,end] = α/2.0*(u[end,end]-u[1,1])\n",
    "    df[end,1:K-1]=α/2.0*(u[end,1:K-1]-u[1,2:K])\n",
    "    mul!(dF_d,D,u)\n",
    "    mul!(dF_c,Φ,df)\n",
    "    return -2.0/h*(dF_d - dF_c)  \n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  1.010782 seconds (2.32 M allocations: 136.950 MiB, 1.99% gc time, 99.95% compilation time)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "retcode: Success\n",
       "Interpolation: 1st order linear\n",
       "t: 2-element Vector{Float64}:\n",
       " 0.0\n",
       " 1.0\n",
       "u: 2-element Vector{Matrix{Float64}}:\n",
       " [0.0 0.5877852522924731 … -0.9510565162951536 -0.5877852522924734; 0.031486749455803056 0.6129671257739798 … -0.9408550126738896 -0.5620204949592386; … ; 0.5620204949592384 0.9408550126738895 … -0.6129671257739798 -0.03148674945580308; 0.5877852522924731 0.9510565162951535 … -0.5877852522924734 -2.4492935982947064e-16]\n",
       " [-0.8414018239068098 -0.36325380133230345 … -0.7738178938452396 -0.9984143611514931; -0.8240522731953155 -0.33365452928529354 … -0.793455206538053 -0.9996505294823455; … ; -0.3923424140010108 0.22325101024290656 … -0.996044018680173 -0.8580630168071015; -0.3631466346382578 0.2537920238242879 … -0.998310871780859 -0.8415061948317838]"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Np = P+1\n",
    "df = zeros(Np, K)\n",
    "dF_d = zeros(Np, K)\n",
    "dF_c = zeros(Np, K)\n",
    "p = (D, Φ, h, df, Np, dF_d, dF_c, K, 0)\n",
    "prob = ODEProblem(rhs_update3, u0, (0.0,1.0),p)\n",
    "@time solve(prob,Tsit5(),save_everystep=false)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "BenchmarkTools.Trial: 7784 samples with 1 evaluation.\n",
       " Range \u001b[90m(\u001b[39m\u001b[36m\u001b[1mmin\u001b[22m\u001b[39m … \u001b[35mmax\u001b[39m\u001b[90m):  \u001b[39m\u001b[36m\u001b[1m407.943 μs\u001b[22m\u001b[39m … \u001b[35m10.693 ms\u001b[39m  \u001b[90m┊\u001b[39m GC \u001b[90m(\u001b[39mmin … max\u001b[90m): \u001b[39m 0.00% … 95.31%\n",
       " Time  \u001b[90m(\u001b[39m\u001b[34m\u001b[1mmedian\u001b[22m\u001b[39m\u001b[90m):     \u001b[39m\u001b[34m\u001b[1m486.758 μs              \u001b[22m\u001b[39m\u001b[90m┊\u001b[39m GC \u001b[90m(\u001b[39mmedian\u001b[90m):    \u001b[39m 0.00%\n",
       " Time  \u001b[90m(\u001b[39m\u001b[32m\u001b[1mmean\u001b[22m\u001b[39m ± \u001b[32mσ\u001b[39m\u001b[90m):   \u001b[39m\u001b[32m\u001b[1m637.130 μs\u001b[22m\u001b[39m ± \u001b[32m 1.285 ms\u001b[39m  \u001b[90m┊\u001b[39m GC \u001b[90m(\u001b[39mmean ± σ\u001b[90m):  \u001b[39m27.12% ± 12.61%\n",
       "\n",
       "  \u001b[34m█\u001b[39m\u001b[32m▂\u001b[39m\u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m▁\n",
       "  \u001b[34m█\u001b[39m\u001b[32m█\u001b[39m\u001b[39m▆\u001b[39m▃\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▆\u001b[39m█\u001b[39m \u001b[39m█\n",
       "  408 μs\u001b[90m        \u001b[39m\u001b[90mHistogram: \u001b[39m\u001b[90m\u001b[1mlog(\u001b[22m\u001b[39m\u001b[90mfrequency\u001b[39m\u001b[90m\u001b[1m)\u001b[22m\u001b[39m\u001b[90m by time\u001b[39m      10.1 ms \u001b[0m\u001b[1m<\u001b[22m\n",
       "\n",
       " Memory estimate\u001b[90m: \u001b[39m\u001b[33m1.78 MiB\u001b[39m, allocs estimate\u001b[90m: \u001b[39m\u001b[33m3362\u001b[39m."
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@benchmark solve(prob,Tsit5(),save_everystep=false)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  0.531198 seconds (1.57 M allocations: 100.317 MiB, 2.82% gc time, 99.79% compilation time)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "retcode: Success\n",
       "Interpolation: 1st order linear\n",
       "t: 2-element Vector{Float64}:\n",
       " 0.0\n",
       " 1.0\n",
       "u: 2-element Vector{Matrix{Float64}}:\n",
       " [0.0 0.5877852522924731 … -0.9510565162951536 -0.5877852522924734; 0.031486749455803056 0.6129671257739798 … -0.9408550126738896 -0.5620204949592386; … ; 0.5620204949592384 0.9408550126738895 … -0.6129671257739798 -0.03148674945580308; 0.5877852522924731 0.9510565162951535 … -0.5877852522924734 -2.4492935982947064e-16]\n",
       " [-0.8414191423804401 -0.3631065856821248 … -0.7739486342866062 -0.9983545188792516; -0.823992142791415 -0.3335934892795848 … -0.7934784812904109 -0.9996510670860366; … ; -0.3922600327893867 0.22332530739193066 … -0.9960361868050959 -0.8580143066961294; -0.3631052455724007 0.25390042410751157 … -0.9983422678399212 -0.8414241644701116]"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@time solve(prob,CVODE_BDF(linear_solver=:GMRES),save_everystep=false)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "BenchmarkTools.Trial: 10000 samples with 1 evaluation.\n",
       " Range \u001b[90m(\u001b[39m\u001b[36m\u001b[1mmin\u001b[22m\u001b[39m … \u001b[35mmax\u001b[39m\u001b[90m):  \u001b[39m\u001b[36m\u001b[1m211.336 μs\u001b[22m\u001b[39m … \u001b[35m 15.714 ms\u001b[39m  \u001b[90m┊\u001b[39m GC \u001b[90m(\u001b[39mmin … max\u001b[90m): \u001b[39m 0.00% … 74.97%\n",
       " Time  \u001b[90m(\u001b[39m\u001b[34m\u001b[1mmedian\u001b[22m\u001b[39m\u001b[90m):     \u001b[39m\u001b[34m\u001b[1m225.879 μs               \u001b[22m\u001b[39m\u001b[90m┊\u001b[39m GC \u001b[90m(\u001b[39mmedian\u001b[90m):    \u001b[39m 0.00%\n",
       " Time  \u001b[90m(\u001b[39m\u001b[32m\u001b[1mmean\u001b[22m\u001b[39m ± \u001b[32mσ\u001b[39m\u001b[90m):   \u001b[39m\u001b[32m\u001b[1m266.805 μs\u001b[22m\u001b[39m ± \u001b[32m709.574 μs\u001b[39m  \u001b[90m┊\u001b[39m GC \u001b[90m(\u001b[39mmean ± σ\u001b[90m):  \u001b[39m10.32% ±  3.81%\n",
       "\n",
       "  \u001b[39m▄\u001b[39m▇\u001b[39m█\u001b[39m▇\u001b[39m▇\u001b[34m█\u001b[39m\u001b[39m▇\u001b[39m▄\u001b[39m▂\u001b[39m▁\u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m▁\u001b[39m \u001b[39m \u001b[39m \u001b[32m \u001b[39m\u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m▂\n",
       "  \u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[34m█\u001b[39m\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m▇\u001b[39m█\u001b[32m█\u001b[39m\u001b[39m█\u001b[39m█\u001b[39m▇\u001b[39m▆\u001b[39m▆\u001b[39m▅\u001b[39m▆\u001b[39m▆\u001b[39m▇\u001b[39m▆\u001b[39m▆\u001b[39m▆\u001b[39m▅\u001b[39m▆\u001b[39m▆\u001b[39m▆\u001b[39m▅\u001b[39m▆\u001b[39m▅\u001b[39m▅\u001b[39m▆\u001b[39m▅\u001b[39m▅\u001b[39m▅\u001b[39m▅\u001b[39m▂\u001b[39m▅\u001b[39m▄\u001b[39m▅\u001b[39m▅\u001b[39m▅\u001b[39m▅\u001b[39m▅\u001b[39m▄\u001b[39m▅\u001b[39m▅\u001b[39m▄\u001b[39m▅\u001b[39m▄\u001b[39m▅\u001b[39m▄\u001b[39m▄\u001b[39m \u001b[39m█\n",
       "  211 μs\u001b[90m        \u001b[39m\u001b[90mHistogram: \u001b[39m\u001b[90m\u001b[1mlog(\u001b[22m\u001b[39m\u001b[90mfrequency\u001b[39m\u001b[90m\u001b[1m)\u001b[22m\u001b[39m\u001b[90m by time\u001b[39m        402 μs \u001b[0m\u001b[1m<\u001b[22m\n",
       "\n",
       " Memory estimate\u001b[90m: \u001b[39m\u001b[33m264.89 KiB\u001b[39m, allocs estimate\u001b[90m: \u001b[39m\u001b[33m1451\u001b[39m."
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@benchmark solve(prob,CVODE_BDF(linear_solver=:GMRES),save_everystep=false)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0005286452943268698"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "uf1 = solve(prob,CVODE_BDF(linear_solver=:GMRES),save_everystep=false)\n",
    "ue = sin.(x .- 1.0)\n",
    "norm(uf1[2]-ue)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0002618762597799035"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "uf2 = solve(prob,Tsit5(),save_everystep=false)\n",
    "ue = sin.(x .- 1.0)\n",
    "norm(uf2[2]-ue)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.4.1",
   "language": "julia",
   "name": "julia-1.4"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.4.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}